Nonadditivity of Fluctuation-Induced Forces in Fluidized Granular Media 
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We investigate the effective long-range interactions between intruder particles immersed in a ran- 
domly driven granular fluid. The effective Casimir-like force between two intruders, induced by 
the fluctuations of the hydrodynamic fields, can change its sign when varying the control parame- 
ters: the volume fraction, the distance between the intruders, and the restitution coefficient. More 
interestingly, by inserting more intruders, we verify that the fluctuation-induced interaction is not 
pairwise additive. The simulation results are qualitatively consistent with the theoretical predictions 
based on mode coupling calculations. These results shed new light on the underlying mechanisms 
of collective behaviors in fluidized granular media. 
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Granular segregation has been extensively investigated 
during the last two decades aimed at revealing the under- 
lying complex dynamics [1, 2]. Besides the scientific in- 
terest, understanding the mechanisms of segregation is of 
essential importance in geophysical [3] and industrial 
processes. The behavior of granular mixtures, when me- 
chanically agitated, depends on a long list of grain, con- 
tainer, and external driving properties [2|,[5|. The control 
parameters can be tuned so that the demixing is initiated, 
reversed, or prevented While the phase behavior 

of these systems is still a matter of debate, the nature of 
particle-particle interactions is known to play a crucial 
role; two extreme limits can be distinguished: (i) the fully 
fluidized regime where particles undergo only binary col- 
lisions, and (ii) the lasting contacts regime where durable 
frictional contacts exist during a considerable part of the 
agitation cycle. While in the latter case the relevant pro- 
cesses are, e.g., reorganization, inertia, and convection 
Uj , some studies reveal the existence of another mecha- 
nism in the fluidized regime: in the presence of intruder 
particles, the hydrodynamic fields are modified especially 
in the inner regions between intruders, leading to effec- 
tive long-range interactions @, [ol, [HI, [l3| . Cattuto et al. 
[l2| found that a pair of intruder particles experience an 
effective force in a driven granular bed, originating from 
the modification of the pressure field fluctuations due to 
the boundary conditions imposed by the intruders. Such 
Casimir-like interactions are expected in thermal n oisy 
environments confined by geometrical constraints |14| . 
Most reports, so far, are about either binary mixtures 
d, Q or one or few intruder particles in a bed of smaller 
ones [1,[iq|5[13|. An important question to address is how 
the collective behavior is influenced by the number and 
arrangement of the intruders. 

In the present Letter, we study the effective interac- 
tions between immobile intruder particles immersed in 
a uniformly agitated granular fluid where all particles 
undergo inelastic binary collisions (Fig. [1]). We show 
that the interaction between a pair of intruders exhibits 



a crossover from attraction to repulsion below a critical 
density, as predicted in 0. We here address the gen- 
eral conditions under which the transition happens, and 
present the phase diagram of the transition. Moreover, 
by comparing the behavior of two and multi intruder sys- 
tems, we find that the fluctuation- induced force is not 
derived from a pair-potential; inserting a new intruder af- 
fects the previously existing interactions in a non-trivial 
way, depending on the relative positions of the intruders. 
Such a feature together with the possible sign change of 
the forces make the multi-body interactions more com- 
plicated and may lead to a variety of collective behaviors 
such as segregation, clustering, or pattern formation. An- 
alytical calculations using the theory of randomly driven 
granular fluids [15] confirm our findings. 

Simulation method — We consider a 2p granular fluid 



similar to the setup described in Refs.[l2|, [iSl, ll6| by 
means of molecular dynamics simulations. We have a 
reference system with two intruders A and B in which 
Lo/r=200, Rjr=l{), and D^^ Jr=?>{) [see Fig. [I^a)]. 
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FIG. 1: (color online) (a) Sketch of the simulation cell, (b) 
^A,B scaled by T^jr vs 0. Comparison is made with the so- 
lution of Eq. ([2]) for the reference system (solid line), as well 
as other values of the control parameters and cx (dashed 
lines), (c) Typical snapshots of dense (gray circles) and dilute 
(red crosses) states with 0=0.66 and 0.24, respectively. 
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FIG. 2: (a) F^^ys D^^ . The solid line is obtained via Eq. ([2|). 
(b) F. 



vs T. The dashed line indicates a linear relation. 
Values of F around lO^^Tg/r reflect the accuracy level of our 
calculation, (c) F^^ vs L. The dashed line corresponds to 
h{L). The inset shows that the deviation of the results from 
h{L) has no systematic dependence on L. 0=0.24 in all cases. 



The reference volume fraction (j)^ is 0.66, and the normal 
restitution coefficient is set to 0.8 for all collisions. Pe- 
riodic boundary conditions are applied in both directions 
of the square-shaped cell to provide a spatially homoge- 
neous state. The system is coupled to an external heat 
bath that uniformly transfers energy into the system; 
the acceleration of each particle is perturbed instanta- 
neously by a random amount which can be considered 
as a Gaussian white noise with zero mean and correla- 
tion {(,ia(t)ijh(t'))=i^5ij5ah^{t — t')^ whcrc a and h denote 
Cartesian components of the vectors, and ^ is the driving 
strength. The rate of the energy gain of a single particle 
averaged over the uncorrelated noise source is dtE=m£^'^ 
[15i] . Each particle also loses energy due to inelastic colli- 
sions at the mean-field rate of dtE={a'^-l)LjT/2 [13, Ell, 
where T is the granular temperature and uj is the collision 
frequency given by the Enskog theory [l9|. Eventually, 
the system reaches a nonequilibrium stationary state by 
balancing the energy input and the dissipation. 

Effective two-body interactions — In the steady state 
we measure the total force exerted by the granular fluid 
on each intruder along the x axis during the time inter- 
val r (r~250 collisions per particle). Due to the ob- 
served large fluctuations, the force is measured for more 
than 10^ consecutive time intervals r. The probabil- 
ity distribution of the data is well fitted by a Gaussian 
[20| with the standard deviation a=0.24:4:T^/r and the 
nonzero mean Fq=— 0.023TQ/r, where is the mean- field 
approximation of the steady-state temperature deduced 
from the Enskog theory [15]. Using a similar analysis 
along the y axis, we obtain zero force within the accu- 
racy of our measurements. E^ can be considered as the 
magnitude of the effective force ^ that the intruder B 
exerts on A, which is attractive in this case. We observe 
that, upon decreasing the volume fraction below a criti- 
cal value (/)^~0.57, the effective interaction ^ becomes 
repulsive [Fig, mb)], in agreement with the prediction 
of Ref. [l2|. However, the transition is controlled not 
only by 0, but also by and a. One expects that. 



far from the transition region, increasing Dj^^ decreases 
the magnitude of ^ and it should eventually vanish at 



D^^=L/2 due to periodic boundary conditions, as con- 
firmed by simulations [Fig.[2]^a)]. In Fig.[2fb), by varying 
the driving strength ^, it is shown that E is proportional 
to the steady state temperature. Moreover, the results 
reveal the impact of dimensionality on the process: E 
increases slightly with L when we vary the system size 
while other parameter values are kept fixed. The simu- 
lation results, shown in Fig. [2jc), can be well fitted by a 
logarithmic growth /i(L)=a ln(L)+6 (dashed line). This 
is contrary to what happens in three dimensional sys- 
tems, where the force is independent of L. 

Triple configurations and nonadditivity — Next we ad- 
dress the interesting case of triple systems, where the 
third intruder is located either on the x or y axis (Fig.[3j). 
By choosing Dj^^/r=30 and 0=0.24, the effective force 
exerted on particle A in the triple configuration 
(A,B,D) is compared to ^ and obtained from the 
binary systems (A,B) and (A,D), respectively. Note that 
the simulation is performed anew for each set of intrud- 
ers. Figure [3fa) shows that is nearly the vector 
sum of i^ g and F^^^. However, a comparison between 
the sets (A,B,C), (A,B), and (A,C) in Figure [3];b) (with 
D^Jr=lQi) reveals that the force is definitely not pairwise 
additive in this case; E^ is even smaller than E^ ^. 

In order to understand the mechanism behind the 
long-range interactions and the transition, we draw at- 
tention to the fact that the fluctuating hydrodynamic 
fields, e.g. density [see Fig. [TJc)], are notably influenced 
by the geometric constraints, resulting in pressure im- 
balance around the intruders and effective interactions 
between them. To establish a quantitative connection 
between the effective force and the hydrodynamic fluc- 
tuations, we first employ mode coupling calculations 
2l| to evaluate the two-body interactions 



12|, [151, |21| to evaluate the two-body interactions. In 
the nonequilibrium steady state, the hydrodynamic fields 
{p{r),T{r)^n{r)) fluctuate around their stationary val- 
ues (ps^Ts^Us). The average pressure fluctuation Pf{r) 
in the presence of the boundary conditions imposed by 
intruders behaves analogously to the Casimir effect, i.e., 
Pf{r) in the hatched region of Fig. [T]^a) differs from that 
ojf the cross-hatched region. Using the Verlet-Levesque 
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FIG. 3: (color online) Gomparison between the binary and 
triple effective forces (scaled by O.OlTo/r) exerted on particle 
A in the presence of particles (a) B and D (b) B and G. (c) e, 
and iA,B+^,c versus ^ position of particle G. 
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FIG. 4: (color online) (a) Schematic phase diagram of the transition in the (0, Dab, cx) space. The surface corresponds to 
3=0. (b-d) 2D profiles of the phase diagram where the solid lines mark the interface position. The dashed lines correspond to 
the interface position for other values of the third control parameter. The color intensity reflects the magnitude of F (scaled by 
O.OITq /r), with blue (red) meaning repulsion (attraction), (e) F on particle A vs n, the length of the chain. The results of the 
mode coupling, simulation, and pairwise summation of two-body forces are shown with solid line, full circles, and open circles, 
respectively. Inset: Mode coupling results (solid line) and their corresponding pairwise summations (dashed line) (L=500r). 



equation of state for a hard disks system p(n, T)=T/(n) 
(with f(n)=n(l + 0^/8)/(l — 0)^, and n the number den- 
sity) [22|, we expand the pressure up to second order 
around (us^Ts)^ and take the statistical average over 
the random noise source: Pj^{r)=f'{n)\^^{Sn{r)ST{r))-\- 
^Tsf" {n)\^^{{5n{r))'^). By Fourier transforming 5n{r) 
and 5T{r) one obtains 



Pf{r)=J 



1 



r{n)\SMk)^-Tsnn)\S^^{k) 



dk, (1) 



where the integral is taken over the k vectors allowed 
at position r by the boundary conditions, and 5'ab(^) is 
the pair structure factor defined as {Sa{k)Sb{—k)) . 
The detailed description of the structure factor calcu- 
lations will be reported elsewhere (see also Ref. [15]). 
Here we denote the integrand of Eq. ^ with g{k^(j)^a) 
and compare Pf{r) for two surface points located on op- 
posite sides of intruder A with the same y coordinates. 
The related k vectors in the x direction are confined to 
Ddy)=D^B-^VlR^-y') and D^Jy)=L- D^,-2^{R^- 
y'^) in the cross-hatched and hatched regions, respec- 
tively; Therefore the pressure difference between these 
two points Ap{y)=p^^^^ — p^^^^^ has y dependence. By 
integrating over y, we arrive at the average pressure dif- 
ference between the gap and outside region: 



rR 

Ap=l dy 

R 



27r/r* p27r/r* 

dk-i dk^ 



^-^-;^^fe,fe,^,0,a)^ (2) 



27r/L 



2R 



To ensure that the hydrodynamic description is valid, 
the integrals are restricted to the long wavelength range 
r*=max(2r, /*) (with P being the mean free path) and 
only small inelasticities are considered. Moreover, is 
always chosen large enough (3i?<I\g) so that the short- 
range depletion forces [23| do not play a role. The ef- 
fective force F^g is calculated via Eq. ^ for different 
values of (j) or and compared to the simulation re- 
sults in Figs. [TJb) and[2fa). The forces are of the same 
order of magnitude as those obtained from the simula- 
tions. The deviations can be attributed to the fact that 



the hydrodynamic fluctuations are correlated in the gap 
and outside regions. The correction due to this effect is 
proportional to d'^g/dcj)'^ which always has a sign opposite 
to that of thus, Eq. (|2]) overestimates the magnitude 
of the force. Using the mode coupHng calculations, it is 
shown in Fig.ITJb) that the transition point is sensitive to 
the choice of and a. The set of control parameters 
for sign switching of the Casimir force crucially depends 
on the physics of the system (see e.g. 0). Figure |4]sum- 
marizes the calculations in a phase diagram in the (0, a, 
^ab) space, which appears to be in remarkable accord 
with the dynamical model. The phase diagram is not 
influenced by the choice of the steady-state temperature 
T5, while the magnitude of the force grows Hnearly with 
Tg as expected for Casimir forces in thermal fluctuating 
media [14]. Regarding the fact, that the leading term 
of g{k^(j)^a) at small k is proportional to 1/k'^ [l5|, one 
also finds from Eq. ([T]) that p^, and therefore the force, 
in the thermodynamic limit behaves as 1/r in 3D while 
diverges logarithmically as ln(L/r) in 2D. 

Multi-body effects — In the triple system of Fig. [31(a), 
loosely speaking, because of the independence of k vec- 
tors in X and y directions one expects that the vector sum 
holds. This is in agreement with the simulation results 
in Fig. [31(a), neglecting the deviation due to hydrody- 
namic correlations. In the configuration of Fig. [31(b), the 
effective force exerted on particle A results from the dif- 
ference between the range of available k modes on its left 
and right sides. In the presence of particle C, the range of 
allowed k modes decreases on the left side due to periodic 
boundary conditions, which causes a lowering of pressure 
difference between both sides of A. The effective force 
is thus smaller, compared to the binary interaction 

3 . The calculated force on particle A using Eq. (|2]) is 
shown with i^ Fig. [31(b); the agreement is satisfac- 

tory. We introduce a measure e=|i^ g+i^ c~^,bc I/I^,b I 
for quantifying the deviation from the case that the in- 
teractions are derived from a pair-potential. Figure [31(c) 
shows how e, obtained from the mode coupling calcula- 
tions, behaves when the position of particle C is varied 



4 



'1'!! 




\ 



FIG. 5: The trajectories of a test intruder subject to the effec- 
tive force field in a fixed square configuration of intruders, ob- 
tained from the mode coupling calculations for 0=0.24 (left), 
and 0.54 (right). Short range interactions are excluded (white 
zones). Lighter colors mean stronger forces. The dashed line 
in the right figure indicates the border of an attractive sub- 
domain, where the trajectories of the test particle end up on 
the fixed intruder. 



on the X axis. Figure HJ^e) indicates that our analytical 
approach also provides reasonable estimates of the effec- 
tive force acting on particle A in a chain configuration, 
while the deviation from pairwise summation of two-body 
forces grows with the number of intruders, n. 

Finally, we investigate more complicated geometries by 
calculating the effective force field of a fixed square struc- 
ture acting on a test intruder particle. To elucidate the 
impact of sign switching on the force pattern we compare 
an intermediate density (around the transition zone) with 
a low density case. Figure [5] shows that the patterns are 
clearly different e.g. in terms of the number of equilibrium 
points. There exist also attractive subdomains of length 
scale ^ around the fixed intruders in the right figure (e.g. 
the domain surrounded by the dashed line). Depending 
on the choice of control parameters, ^ ranges between 
(entirely repulsive patterns at low densities, as shown in 
the left figure) and L (entirely attractive patterns at high 
densities, not shown) leading to different types of collec- 
tive behavior such as segregation, clustering, and pattern 
formation. Putting aside the pairwise summations, we 
compare the mode coupling predictions with simulation 
results at four selected points in Fig [51 The predicted 



TABLE L Comparison between the theoretical (T) and sim- 
ulation (S) results for the selected points in Figs. [5jleft) and 
[5{right). Same force scales as Fig.[3l The maximum error bars 
for the simulated F and are zbO.l and ±7°, respectively. 



\F\{A) \F\{B) \F\{C) \F\{D) 0{A) 0(B) 



left-S 1.0 2.2 3.4 3.7 13° 92° 51° 44° 

left-T 1.17 2.52 3.97 4.32 7° 90° 60° 45° 

right-S 0.2 0.2 0.2 0.0 19° 88° -32° - 

right-T 0.24 0.26 0.25 0.02 28° 90° -24° 225° 



forces are of the same order of magnitude as the simula- 
tion results (see table H]), however, the errors of the force 
size \F\ and its direction 9 reach up to 30% and 15°, re- 
spectively. The differences are smaller at large distances 
and also low densities. Indeed, the hydrodynamic corre- 
lations become more important in multi-body cases and, 
hence, one should include triplet or higher order struc- 
ture factors [m in mode coupling calculations to properly 
take the multi-body effects into account. 

In conclusion, we focus on the problem of long-range 
fluctuation-induced forces between intruder particles im- 
mersed in an agitated fluid bed. The sign of the force 
can be reversed by tuning the control parameters. Fur- 
thermore, the multi-body interactions do not follow from 
two-body force descriptions. A newly inserted intruder, 
depending on its position, may affect the pressure balance 
around the other intruders. This suggests that the effec- 
tive force is not derived from a pair-potential in agree- 
ment with our simulation results. Our findings represent 
a step forward in understanding the origin of collective 
behaviors in fluidized granular mixtures. 
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